Figure 2

The Figure 2 is the showcase for geomeTriD package to present multiple genomic signals along with 3D models for a microcompartment/loop.

Load Libraries

library(geomeTriD)
library(geomeTriD.documentation)
library(GenomicRanges)
library(TxDb.Mmusculus.UCSC.mm39.knownGene)
library(org.Mm.eg.db)
library(trackViewer)
library(RColorBrewer)

Prepare the annotation and genomic data

The Region Capture Micro-C data were downloaded from GSE207225. The genomic signals of ChIP-seq were downloaded from GSE178982 and remapped to mm39 genome.

## import Genomic interaction data 
geo_acc <- c("DMSO"='GSM6281851',
             'IAA'='GSM6281852')
url <- 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn';
urls <- mapply(geo_acc, names(geo_acc), FUN=function(.f, .cond){
  file.path(url, .f, 'suppl', 
            paste0(.f, '_RCMC_BR1_merged_allCap_', .cond,
                   '_mm39.merged.50.mcool'),
            fsep='/')
})
## import the genomic interactions
range_chr8_wid <- GRanges('chr8:84000000-92000000')
gi <- importGInteractionsFromUrl(urls=urls, resolution=500, range=range_chr8_wid,
                                 format='cool', normalization='balanced')
## adding rname 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn/GSM6281851/suppl/GSM6281851_RCMC_BR1_merged_allCap_DMSO_mm39.merged.50.mcool'
## adding rname 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn/GSM6281852/suppl/GSM6281852_RCMC_BR1_merged_allCap_IAA_mm39.merged.50.mcool'
## create 3d model by mdsPlot function
range_chr8 <-  GRanges('chr8:85550000-85800000')
model3d <- lapply(gi, mdsPlot, range = range_chr8, k=3, render = 'granges')
## Warning in checkGI(gi, fixedBin = TRUE): There are NA values in the gi score.
## It will be removed.
## Some region are missing from the input gi.
## initial  value 17.000699 
## iter   5 value 14.159797
## iter  10 value 13.013113
## iter  15 value 12.575326
## final  value 12.445963 
## converged
## Warning in checkGI(gi, fixedBin = TRUE): There are NA values in the gi score.
## It will be removed.
## Some region are missing from the input gi.
## initial  value 17.114933 
## iter   5 value 13.719667
## iter  10 value 12.444197
## iter  15 value 12.102132
## final  value 11.986303 
## converged
### create gene annotations
#### get all genes
feature.gr <- genes(TxDb.Mmusculus.UCSC.mm39.knownGene)
##   2 genes were dropped because they have exons located on both strands of the
##   same reference sequence or on more than one reference sequence, so cannot be
##   represented by a single genomic range.
##   Use 'single.strand.genes.only=FALSE' to get all the genes in a GRangesList
##   object, or use suppressMessages() to suppress this message.
#### subset the data by target viewer region
feature.gr <- subsetByOverlaps(feature.gr, range_chr8)
#### assign symbols for each gene
symbols <- mget(feature.gr$gene_id, org.Mm.egSYMBOL, ifnotfound = NA)
feature.gr$label[lengths(symbols) == 1] <- 
  unlist(symbols[lengths(symbols) == 1])
#### assign colors for each gene
feature.gr$col <- sample(1:7, length(feature.gr), replace = TRUE)

### plot for the detailed region ii to show the open or close loop.
subregion <- GRanges('chr8:85690000-85730000')
model3d.sub <- lapply(model3d, subsetByOverlaps, ranges=subregion)
model3d.sub[['IAA']] <- alignCoor(model3d.sub[['IAA']], model3d.sub[['DMSO']])

### import cohesion ctcf, and yy1 signals, realigned to mm39 for GSE178982
pf <- system.file('extdata', package = 'geomeTriD.documentation')
bws_files <- dir(file.path(pf, 'GSE178982'), '.bw', full.names = TRUE)
(names(bws_files) <- sub('^.*?_(IAA|UT)_(.*?).CPM.*$', '\\2', bws_files))
##  [1] "CTCF"  "RAD21" "SMC1A" "SMC3"  "YY1"   "CTCF"  "RAD21" "SMC1A" "SMC3" 
## [10] "YY1"
bw_UT <- bws_files[grepl('_UT_', bws_files)]
bw_IAA <- bws_files[grepl('_IAA_', bws_files)]
signals_UT <- lapply(bw_UT, importScore, format='BigWig', ranges=subregion)
signals_IAA <- lapply(bw_IAA[names(bw_UT)], importScore, format='BigWig', ranges=subregion)
colorSets <- c(CTCF="cyan",YY1="yellow",
               RAD21="red", SMC1A="green", SMC3="blue")
for(i in seq_along(signals_UT)){
  setTrackStyleParam(signals_UT[[i]], "color", c("gray30", colorSets[names(signals_UT)[i]]))
  setTrackStyleParam(signals_IAA[[i]], "color", c("gray30", colorSets[names(signals_UT)[i]]))
}
## set maximal lwd for UT and IAA samples according their max
genomicScoreRanges <- lapply(signals_UT, function(.ele) range(.ele$dat$score))

Plot the data by geomeTriD

### create the structure
dmso <- view3dStructure(model3d.sub[['DMSO']], feature.gr = feature.gr,
                        genomicSigs=signals_UT,
                        reverseGenomicSigs = FALSE,
                        genomicScoreRange = genomicScoreRanges,
                        lwd.maxGenomicSigs = 20,
                        k = 3, renderer = 'none')
## Feature type is missing. Set as gene.
iaa <- view3dStructure(model3d.sub[['IAA']], feature.gr = feature.gr,
                       genomicSigs=signals_IAA,
                       reverseGenomicSigs = FALSE,
                       genomicScoreRange = genomicScoreRanges,
                        lwd.maxGenomicSigs = 20,
                       k=3, renderer = 'none')
## Feature type is missing. Set as gene.
iaa <- lapply(iaa, function(.ele){
  .ele$side <- 'right'
  .ele
})
threeJsViewer(dmso, iaa, title = c('DMSO control', 'IAA 3h'))
#widget <- threeJsViewer(dmso, iaa, title = c('DMSO control', 'IAA 3h'), background = 'white')
#tempfile <- 'Fig2.part1.html'#tempfile(fileext = '.html', pattern = 'RCMC_BR1_IAA_vs_DMSO_Klf1_II.3jsViewer.')
#htmlwidgets::saveWidget(widget, file=tempfile)
#utils::browseURL(tempfile)

Plot with interaction signals

The following plot clearly highlights the top 10 interaction events for DMSO and IAA samples. These interactions shift from the long-distance regions of Hook2 to more localized regions.

## extract backbone coordinates, which will be used as the bone for RCMC data
xyz_dmso <- extractBackbonePositions(dmso)
xyz_iaa <- extractBackbonePositions(iaa)

DMSO <- gi[['DMSO']]
IAA <- gi[['IAA']]

DMSO.subset_II <- subsetByOverlaps(DMSO, subregion)
IAA.subset_II <- subsetByOverlaps(IAA, subregion)
DMSO.subset_II <- DMSO.subset_II[distance(first(DMSO.subset_II), second(DMSO.subset_II))>5000]
IAA.subset_II <- IAA.subset_II[distance(first(IAA.subset_II), second(IAA.subset_II))>5000]
hic_dmso_II <- create3dGenomicSignals(
  DMSO.subset_II, 
  xyz_dmso,
  name='dmsoII', # name prefix for the geometry
  tag='dmsoII', # name for the layer in the scene
  color = c('white', brewer.pal(9, 'YlOrRd')),
  topN=10, # only plot the top 10 events ordered by the scores.
  lwd.maxGenomicSigs = 3,
  alpha=0.5
)
hic_iaa_II <- create3dGenomicSignals(
  IAA.subset_II, 
  xyz_iaa,
  name='iaaII', # name prefix for the geometry
  tag='iaaII', # name for the layer in the scene
  color = c('white', brewer.pal(9, 'YlOrRd')),
  topN=10, # only plot the top 10 events ordered by the scores.
  lwd.maxGenomicSigs = 3,
  alpha=0.5
)
hic_iaa_II <- lapply(hic_iaa_II, function(.ele){
  .ele$side <- 'right'
  .ele
})
threeJsViewer(dmso, iaa, hic_dmso_II, hic_iaa_II, title = c('DMSO', 'IAA 3h'), background = c('gray10', 'gray20', 'gray20', 'gray10'))
#widget <- threeJsViewer(dmso, iaa, hic_dmso_II, hic_iaa_II, title = c('DMSO', 'IAA 3h'), background = c('gray10', 'gray20', 'gray20', 'gray10'))
#tempfile <- 'Fig2.part2.html'#tempfile(fileext = '.html', pattern = 'RCMC_BR1_IAA_vs_DMSO_Klf1_II.3jsViewer.')
#htmlwidgets::saveWidget(widget, file=tempfile)
#utils::browseURL(tempfile)

SessionInfo

## R Under development (unstable) (2025-01-08 r87545)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-3                       
##  [2] trackViewer_1.43.6                       
##  [3] org.Mm.eg.db_3.20.0                      
##  [4] TxDb.Mmusculus.UCSC.mm39.knownGene_3.20.0
##  [5] GenomicFeatures_1.59.1                   
##  [6] AnnotationDbi_1.69.0                     
##  [7] Biobase_2.67.0                           
##  [8] GenomicRanges_1.59.1                     
##  [9] GenomeInfoDb_1.43.2                      
## [10] IRanges_2.41.2                           
## [11] S4Vectors_0.45.2                         
## [12] BiocGenerics_0.53.3                      
## [13] generics_0.1.3                           
## [14] geomeTriD.documentation_0.0.1            
## [15] geomeTriD_1.1.6                          
## 
## loaded via a namespace (and not attached):
##   [1] strawr_0.0.92               rstudioapi_0.17.1          
##   [3] jsonlite_1.8.9              magrittr_2.0.3             
##   [5] rmarkdown_2.29              fs_1.6.5                   
##   [7] BiocIO_1.17.1               zlibbioc_1.53.0            
##   [9] ragg_1.3.3                  vctrs_0.6.5                
##  [11] memoise_2.0.1               Rsamtools_2.23.1           
##  [13] RCurl_1.98-1.16             base64enc_0.1-3            
##  [15] htmltools_0.5.8.1           S4Arrays_1.7.1             
##  [17] progress_1.2.3              plotrix_3.8-4              
##  [19] curl_6.1.0                  Rhdf5lib_1.29.0            
##  [21] rhdf5_2.51.2                SparseArray_1.7.3          
##  [23] Formula_1.2-5               sass_0.4.9                 
##  [25] bslib_0.8.0                 htmlwidgets_1.6.4          
##  [27] desc_1.4.3                  Gviz_1.51.0                
##  [29] httr2_1.0.7                 cachem_1.1.0               
##  [31] GenomicAlignments_1.43.0    igraph_2.1.3               
##  [33] lifecycle_1.0.4             pkgconfig_2.0.3            
##  [35] Matrix_1.7-1                R6_2.5.1                   
##  [37] fastmap_1.2.0               GenomeInfoDbData_1.2.13    
##  [39] MatrixGenerics_1.19.1       digest_0.6.37              
##  [41] colorspace_2.1-1            textshaping_0.4.1          
##  [43] Hmisc_5.2-2                 RSQLite_2.3.9              
##  [45] filelock_1.0.3              httr_1.4.7                 
##  [47] abind_1.4-8                 compiler_4.5.0             
##  [49] withr_3.0.2                 bit64_4.6.0-1              
##  [51] htmlTable_2.4.3             backports_1.5.0            
##  [53] BiocParallel_1.41.0         DBI_1.2.3                  
##  [55] biomaRt_2.63.0              MASS_7.3-64                
##  [57] rappdirs_0.3.3              DelayedArray_0.33.3        
##  [59] rjson_0.2.23                tools_4.5.0                
##  [61] foreign_0.8-87              nnet_7.3-20                
##  [63] glue_1.8.0                  restfulr_0.0.15            
##  [65] InteractionSet_1.35.0       rhdf5filters_1.19.0        
##  [67] checkmate_2.3.2             cluster_2.1.8              
##  [69] gtable_0.3.6                BSgenome_1.75.0            
##  [71] ensembldb_2.31.0            data.table_1.16.4          
##  [73] hms_1.1.3                   xml2_1.3.6                 
##  [75] XVector_0.47.2              pillar_1.10.1              
##  [77] stringr_1.5.1               dplyr_1.1.4                
##  [79] BiocFileCache_2.15.0        lattice_0.22-6             
##  [81] deldir_2.0-4                rtracklayer_1.67.0         
##  [83] bit_4.5.0.1                 biovizBase_1.55.0          
##  [85] tidyselect_1.2.1            Biostrings_2.75.3          
##  [87] knitr_1.49                  gridExtra_2.3              
##  [89] ProtGenerics_1.39.2         SummarizedExperiment_1.37.0
##  [91] xfun_0.50                   matrixStats_1.5.0          
##  [93] stringi_1.8.4               UCSC.utils_1.3.1           
##  [95] lazyeval_0.2.2              yaml_2.3.10                
##  [97] evaluate_1.0.3              codetools_0.2-20           
##  [99] interp_1.1-6                tibble_3.2.1               
## [101] cli_3.6.3                   rpart_4.1.24               
## [103] systemfonts_1.2.0           munsell_0.5.1              
## [105] jquerylib_0.1.4             dichromat_2.0-0.1          
## [107] Rcpp_1.0.14                 grImport_0.9-7             
## [109] dbplyr_2.5.0                png_0.1-8                  
## [111] XML_3.99-0.18               parallel_4.5.0             
## [113] pkgdown_2.1.1               rgl_1.3.16                 
## [115] ggplot2_3.5.1               blob_1.2.4                 
## [117] prettyunits_1.2.0           jpeg_0.1-10                
## [119] latticeExtra_0.6-30         AnnotationFilter_1.31.0    
## [121] bitops_1.0-9                txdbmaker_1.3.1            
## [123] VariantAnnotation_1.53.1    scales_1.3.0               
## [125] purrr_1.0.2                 crayon_1.5.3               
## [127] rlang_1.1.5                 KEGGREST_1.47.0